Land Cover 2010 to 2020#

From 2010 to 2020, the Denver MSA has maintained control of urban sprawl. While the population increased by 16%, from 2.54 million to 2.95 million (block-level population from the 2010 and 2020 decennial census datasets), the urban area has expanded by only 3.5%, from 1,200 to 1,240 square miles, according to data from MRLC.

Hide code cell content
import numpy as np
import pandas as pd
import geopandas as gpd
import rasterio as rio
from shapely.geometry import Polygon, Point
crs = "EPSG:2232"
%run assets/colors.py

We used land cover data from 2008 and 2019 as proxies of the land cover in 2010 and 2020, as these are the best available data from MRLC. Land cover data rasters with 4000*4000 sqft cell sizes. Our study identified 5 types of land cover types:

  • Urban.

  • Forests.

  • Farms.

  • Wetlands.

  • Water.

Hide code cell content
# Read land cover rasters
# The file paths are relative to the notebook source file's location in the repo's main branch
land_cover_2009 = rio.open("../data/land-cover/2009.tif")
land_cover_2019 = rio.open("../data/land-cover/2019.tif")

def raster_to_shape(raster, crs, to_point=False, value_column="value"):
    data_array = raster.read(1)
    transform = raster.transform
    shapes_data = []
    height, width = data_array.shape

    if to_point:
        for row in range(height):
            for col in range(width):
                value = data_array[row, col]
                if np.isnan(value):
                    continue
                x, y = transform * (col + 0.5, row + 0.5)
                # Create a polygon for the cell
                point = Point(x, y)
                # Add the polygon and its value to the list
                shapes_data.append({"geometry": point, value_column: value})
        return gpd.GeoDataFrame(shapes_data, crs=crs)

    for row in range(height):
        for col in range(width):
            value = data_array[row, col]
            if np.isnan(value):
                continue
            x_min, y_max = transform * (col, row)
            x_max, y_min = transform * (col + 1, row + 1)
            # Create a polygon for the cell
            polygon = Polygon(
                [(x_min, y_min), (x_min, y_max), (x_max, y_max), (x_max, y_min)]
            )
            # Add the polygon and its value to the list
            shapes_data.append({"geometry": polygon, value_column: value})
    return gpd.GeoDataFrame(shapes_data, crs=crs)

# Get a fishnet of info from the land cover rasters
fishnet = (
    raster_to_shape(
        land_cover_2009, crs, to_point=False, value_column="land_cover_2009"
    )
    .query("land_cover_2009 != 0")
    .copy()
)
points_2019 = (
    raster_to_shape(land_cover_2019, crs, to_point=True, value_column="land_cover_2019")
    .query("land_cover_2019 != 0")
    .copy()
)
# Join points to fishnet, no index right
fishnet = gpd.sjoin(fishnet, points_2019, how="left", predicate="contains").drop(
    columns=["index_right"]
)

# Classify land cover types
developed = [21, 22, 23, 24]
land_covers = {
    "developed": [21, 22, 23, 24],
    "forest": [41, 42, 43],
    "farm": [81, 82],
    "wetland": [90, 95],
    "water": [11],
}
for year in [2009, 2019]:
    # First update whether developed bool
    fishnet[f"developed_{year}"] = fishnet[f"land_cover_{year}"].isin(developed)

    # Then update land cover
    fishnet[f"land_cover_type_{year}"] = "other"
    for land_cover, codes in land_covers.items():
        fishnet.loc[
            fishnet[f"land_cover_{year}"].isin(codes), f"land_cover_type_{year}"
        ] = land_cover

    fishnet.drop(columns=[f"land_cover_{year}"], inplace=True)
Hide code cell source
import altair as alt
alt.data_transformers.disable_max_rows()

def altair_fishnet(fishnet, column, color_dict, legend_title, title):
    chart = alt.Chart(
        fishnet.to_crs(4326)[["geometry", column]]
    ).mark_geoshape().encode(
        color=alt.Color(
            f"{column}:N",
            title=legend_title,
            scale=alt.Scale(
                domain=list(color_dict.keys()), range=list(color_dict.values())
            ),
        )
    ).properties(
        width=400, height=330, title=title
    )
    return chart

# Assets is in the source folder in the repo
from assets.colors import (
    palette_primary,
    palette_backup2_dimmed,
    palette_backup1,
    palette_green,
    palette_green_faded,
    palette_water,
    palette_hero,
)
color_dict = {
    "developed": palette_hero,
    "forest": palette_backup2_dimmed,
    "farm": palette_green,
    "wetland": palette_green_faded,
    "water": palette_water,
    "other": "#dddddd",
}

altair_fishnet(
    fishnet=fishnet,
    column="land_cover_type_2009",
    color_dict=color_dict,
    legend_title="Land Cover Type",
    title="Land Cover in 2008",
)
Hide code cell source
altair_fishnet(
    fishnet=fishnet,
    column="land_cover_type_2019",
    color_dict=color_dict,
    legend_title="Land Cover Type",
    title="Land Cover in 2019",
)
Hide code cell source
fishnet["change"] = "Undeveloped"
fishnet.loc[fishnet["developed_2009"], "change"] = "Developed by 2010"
fishnet.loc[
    (fishnet["developed_2019"] & (~fishnet["developed_2009"])), "change"
] = "Developed 2010-2020"

color_dict = {
    "Developed by 2010": palette_hero,
    "Developed 2010-2020": palette_primary,
    "Undeveloped": "#dddddd",
}

altair_fishnet(
    fishnet=fishnet,
    column="change",
    color_dict=color_dict,
    legend_title="Development Change",
    title="Development Change 2010-2020",
)

The above plots show that the Denver MSA’s urban growth was modest during 2010-2020. Most new developments were sporadically located on the fringes of existing development. We can also see a patter that developments tend to be located near some key corridors to the southeast of the City of Denver. We will further explore these patterns in the next sections.

Hide code cell content
fishnet.to_file("../data/fishnets/with_land_cover.geojson", driver="GeoJSON")